#include <stdio.h>
#include <stdlib.h>
#include <gmp.h>
#include <time.h>
#include <pthread.h>

#define ITER 35 //Numero de iteracoes que serao executadas

//Declaracao das variaveis utilizadas durante as iteracoes. Declaradas publicamente para que as threads tenham acesso.
mpf_t a[ITER]; 
mpf_t b[ITER];
mpf_t p[ITER];
pthread_mutex_t fA[ITER];
pthread_mutex_t fB[ITER];

//Funcoes que calculam os valores de cada variavel a cada iteracao
void *f_calculaA(void *empty) 
{ 
	long i;
	mpf_t aux;
	mpf_init2(aux, 40000000);
	for (i = 1 ; i < ITER ; i++)
	{
		mpf_sqrt(aux, a[i-1]);
		mpf_ui_div(a[i], 1, aux);
		mpf_add(a[i], a[i], aux);
		mpf_div_ui(a[i], a[i], 2); //an+1
		pthread_mutex_unlock(&fA[i]);
	}
}

void *f_calculaB(void *empty)
{
	mpf_t aux;
	mpf_init2(aux, 40000000);
	long i;
	for (i = 1 ; i < ITER ; i++)
	{
		pthread_mutex_lock(&fA[i-1]);
		pthread_mutex_unlock(&fA[i-1]);
		mpf_sqrt(aux, a[i-1]);
		mpf_add_ui(b[i], b[i-1], 1);
		mpf_mul(b[i], b[i], aux);
		mpf_add(aux, a[i-1], b[i-1]);
		mpf_div(b[i], b[i], aux); //bn+1
		pthread_mutex_unlock(&fB[i]);
	}
}

void *f_calculaP(void *empty)
{
	mpf_t aux;
	mpf_init2(aux, 40000000);
	long i;
	for (i = 1 ; i < ITER ; i++)
	{
		pthread_mutex_lock(&fA[i]);
		pthread_mutex_unlock(&fA[i]);
		pthread_mutex_lock(&fB[i]);
		pthread_mutex_unlock(&fB[i]);
		mpf_add_ui(p[i], a[i], 1);
		mpf_mul(p[i], p[i], p[i-1]);
		mpf_mul(p[i], b[i], p[i]);
		mpf_add_ui(aux, b[i], 1);
		mpf_div(p[i], p[i], aux);
	}
}		

int main (int argc, char *argv[]) 
{
	//Variaveis usadas para calcular tempo levado pelo algoritmo
	time_t begin, end;
	double time_spent;
	
	//Declaracao de cada uma das threads que serao executadas
	pthread_t calculaA;
	pthread_t calculaB;
	pthread_t calculaP;
	
	long i;
	mpf_t aux;
	
	
	//Abertura do arquivo onde sera armazenado o valor final calculado
	FILE *piFile = fopen("out_borweinPar.txt","w+");
	if(piFile == NULL) 
	{
		return EXIT_FAILURE;
	}
	
	//Setar precisão (usar todas as 10 milhões de casas) e inicializar variaveis do metodo
	for (i = 1 ; i < ITER ; i++)
	{
		mpf_init2(a[i], 40000000);
		mpf_init2(b[i], 40000000);
		mpf_init2(p[i], 40000000);
		pthread_mutex_init(&fA[i],NULL);
		pthread_mutex_lock(&fA[i]);
		pthread_mutex_init(&fB[i],NULL);
		pthread_mutex_lock(&fB[i]);
	}
	
	mpf_init2(aux, 40000000);
	mpf_init2(a[0], 40000000);
	mpf_init2(b[0], 40000000);
	mpf_init2(p[0], 40000000);
	mpf_sqrt_ui(aux, 2); //a0 = 2^(1/2)
	mpf_set(a[0], aux);
	mpf_set_ui(b[0], 0); //b0 = 0
	mpf_add_ui(p[0], aux, 2); //p0 = 2 + 2^(1/2) 
	
	printf("Starting...\n");
	
	//Inicializa o cronometro
	time(&begin);
	
	pthread_create(&calculaA,NULL,f_calculaA,NULL);
	pthread_create(&calculaB,NULL,f_calculaB,NULL);
	pthread_create(&calculaP,NULL,f_calculaP,NULL);
	
	pthread_join(calculaA,NULL);
	pthread_join(calculaB,NULL);
	pthread_join(calculaP,NULL);
	

	//Finaliza o cronometro e mostra na tela o tempo levado pelo programa
	time(&end);
	time_spent = difftime(end, begin);
	printf("\nExecution: %lf s\nNumber of Iterations: %ld\n",time_spent,i);

	//Escreve no arquivo o valor de pi
	mpf_out_str(piFile,10,250000,p[ITER-1]);
	
	//Fecha o arquivo e limpa as variaveis do algoritmo
	fclose(piFile);
	//mpf_clears(aN,bN,pN,aNP1,bNP1,pNP1,NULL);
}
